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Research  Objectives 


As  stated  in  the  research  proposal,  the  objective  of  this  program  was  to  perform  basic 
research  on  the  application  of  estimation-theoretic  processing  methods  for  unconventional 
imaging  systems.  Imaging  applications  involving  both  incoherent  (passive)  and  coherent 
(active)  illumination  were  to  be  investigated,  and  the  four  significant  thrusts  of  the  program 
were  to  be: 

1.  the  development  of  optimal  processing  methods  for  identifying  and  correcting  wave- 
front  errors  that  arise  due  to  a  non-homogeneous  propagation  medium  or  sensor- 
platform  instabilities; 

2.  the  utilization  of  statistical  priors  and  penalties  in  the  image-estimation  and  system- 
identification  procedures; 

3.  the  utilization  of  novel  and/or  auxiliary  measurements  such  as  those  supplied  by  a 
wave-front  sensor  or  pupil  mask;  and 

4.  the  development  of  improved  methods  for  image  formation  from  sparse  pupil-plane 
arrays,  and  an  investigation  into  the  optimal  number  and  placement  of  array  elements 
for  coherent  active-illumination  systems. 

Because  funding  for  this  project  was  reduced  from  three  to  two  years,  most  of  the  efforts 
were  restricted  to  the  first  three  objectives.  Collaborative  efforts  were  conducted  with  Air 
Force  personnel  at  the  Air  Force  Research  Laboratory  in  Albuquerque,  NM  and  the  Air 
Force  Maui  Optical  Station  in  Maui,  HI. 
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Research  Methods  and  Results 


This  section  documents  the  research  methods  used  and  results  obtained  for  an  investigation 
of  estimation-theoretic  methods  for  unconventional  imaging.  The  methods  and  results  are 
summarized  for  two  areas: 

1.  Estimation-theoretic  bounds  on  the  estimation  accuracy  for  multi-channel  phase  re¬ 
trieval;  and 

2.  Dual-channel  signal  recovery  from  auto-  and  cross-correlation  data. 

Lower  Bounds  on  Estimation  Accuracy  for  Multi-Channel  Phase 
Retrieval 

Wavefield  aberrations  induced  by  optical  defects  in  monolithic  lenses  or  mirrors,  alignment  er¬ 
rors  in  segmented-aperture  systems,  or  wave  propagation  through  an  inhomogeneous  medium 
can  severely  limit  the  resolution  of  an  optical  imaging-system  and  result  in  significantly  dis¬ 
torted  imagery.  Adaptive  compensation  or  post-detection  processing  can  be  used  to  correct 
for  the  aberrations  and  produce  undistorted  imagery,  but  these  procedures  require  accurate 
estimation  of  the  wavefield  aberrations  -  a  phase-retrieval  problem  must  be  solved. 

Hartmann  sensors  are  commonly  used  to  estimate  wavefront  or  optical  aberrations  in  an 
imaging  system.  These  sensors  operate  on  the  principle  that  the  wavefront  aberrations  are 
well-approximated  by  linear  tilts  over  small  spatial  regions  -  a  Hartmann  sensor  segments 
the  wavefield  into  small  spatial  segments  with  a  lenslet  array,  and,  in  principle,  the  spots 
formed  by  the  array  translate  according  to  the  local  wavefront  slope  [1].  The  local  slope 
estimates  are  then  processed  to  estimate  the  wavefront  (or  aberrations).  The  point-spread 
function  from  a  conventional  image  can  also  be  used  to  estimate  the  wavefront  or  aberrations, 
but  this  is  done  less  frequently.  This  method  was  used  out  of  necessity,  however,  to  analyze 
the  aberrations  of  the  Hubble  Space  Telescope’s  infamous  primary  mirror  [2].  Another 
method  that  can  be  applied  to  the  wavefront  or  optical-aberration  estimation  problem  is 
phase  diversity,  whereby  part  of  the  light  from  a  conventional  image  is  split  into  a  diversity 
channel  where  an  additional  aberration  is  introduced.  This  aberration  is  frequently  induced 
by  collecting  the  diversity  data  out  of  the  focal  plane  [3] . 

Performance  comparisons  for  the  wavefront-  and  aberration-estimation  systems  men¬ 
tioned  above  are  typically  linked  to  specific  processing  schemes.  In  contrast  with  this 
approach,  our  goal  in  this  project  is  to  provide  the  fundamental  limits  on  wavefront-  or 
aberration-estimation  accuracy  that  are  imposed  by  the  use  of  various  sensing  modalities. 
To  accomplish  this,  the  Fisher  information  and  Cramer-Rao  bound  are  used  to  quantify  the 
information  content  of  various  modalities,  and  simulation  results  are  presented  to  compare 
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optimal  estimation  performance  to  the  bounds.  A  surprising  result  of  the  analysis  shows  that 
data  recorded  by  a  conventional  imaging  system  contains  more  information  for  parameter 
identification  than  do  those  recorded  with  a  Hartmann  sensor. 


Mathematical  Models 


For  quasi-monochromatic  wave  propagation  through  a  relatively  homogeneous  medium,  the 
diffraction-limited  incoherent  point-spread  function  for  a  conventional  imaging  system  or 
Hartmann  wavefront  sensor  is  well-modeled  as: 


Hv) 


/  T(u)s(y,u)du 
Jv 


(1) 


where  y  and  u  are  two-dimensional  position  variables  in  the  sensor  and  pupil  planes,  respec¬ 
tively,  T(-)  is  the  pupil  transmittance  function,  V  is  the  pupil  region,  and  s(-,  •)  is  the  system 
transfer  function  for  propagation  from  the  pupil  plane  to  the  sensor  plane.  For  the  case  of 
an  in-focus  imaging  system  in  the  Fresnel  approximation,  the  pupil  transmittance  function 
is  modeled  as 

T(u)  =  A(u) exp  ,  (2) 

and  the  system  propagation  function  is  modeled  as 


s(y,u)  = 


exp  ( j  ^fdi)  f.n\y-u\ 


jXdt 


exp 


Xdi 


(3) 


where  A  is  the  nominal  wavelength  of  the  detected  radiation,  /  is  the  system  focal  length, 
A(-)  is  the  system  aperture  or  pupil  function,  and  di  is  the  distance  from  the  pupil  plane  to 
sensor  plane.  For  a  P-element  Hartmann  type  lenslet  array,  the  pupil  transmittance  function 
can  be  modeled  as 

T(u)  =  ^Al(u-up)exv^-j-^\u-up\^J  ,  (4) 

where  up  denotes  the  location  of  the  pth  lenslet  in  the  array,  Ai(-)  is  the  lenslet  aperture 
function,  and  /*  is  the  lenslet  focal  length.  In  this  situation  the  Fresnel  approximation  for 
propagation  from  the  lenslet  plane  to  the  sensor  plane  may  not  be  adequate,  and  the  more 
general  Huygens-Fresnel  propagation  function: 


s(y,u) 


1  exp(j^|?/--ul 

\J\y  -  u?  + 


l  +  ti) 

dj 
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should  be  used,  where  di  is  the  distance  from  the  lenslet  plane  to  the  sensor  plane. 


(5) 
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For  situations  involving  optical  aberrations  in  monolithic  lenses  or  mirrors,  alignment 
errors  in  segmented  systems,  or  wave  propagation  through  an  inhomogeneous  medium,  the 
system  point-spread  function  is  often  modeled  with  a  single  phase-screen  aberration: 


h(y;9)  =  J^T(u)exp[jO(u)\s(y,u)du 


(6) 


where  O(-)  denotes  the  phase-error  (in  radians)  caused  by  the  system  aberration,  and  T(-) 
is  the  abberation-free  pupil  transmittance  function.  The  strength  of  an  aberration  is  often 
quantified  through  the  average  of  its  square  over  the  pupil  region  V: 


e 
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(7) 


but  care  must  be  used  when  applying  this  measure  because  constant  offsets  (piston)  and 
linear  terms  (tilt)  will  increase  e2,  but  will  not  affect  the  quality  of  a  system’s  imagery. 

For  systems  that  divide  the  incoming  light  into  multiple  channels  -  such  as  phase  diversity 
or  the  simultaneous  collection  of  a  conventional  image  and  a  Hartmann  wavefront-sensor 
image  -  the  aberrated  point-spread  functions  for  each  channel  can  be  modeled  as: 


hk{y\ 0)=Oik  Tk(u )  exp  [jd(u)}  sk(y,  u)du 


k  =  l,2,...,K, 


(8) 


where  the  channel  point-spread  functions  are  indexed  by  k  and  each  channel  is  characterized 
by  its  pupil  transmittance  function  Tk(-)  and  propagation  kernel  s*:(-,-).  The  scale  factor 
ak  accounts  for  both  the  exposure  time  and  the  percentage  of  received  light  that  is  diverted 
to  the  kth  channel.  For  situations  of  interest  here,  the  transmittance  functions  for  each 
channel  Tk(-)  can  be  normalized  so  that  hk{y\6 )  integrates  to  ak.  Then  ak  is  a  parameter 
that  determines  the  expected  number  of  photocounts  (light  level)  in  the  kth  channel. 

In  many  applications  and  analyses  it  is  common  to  expand  the  wavefront  error  in  terms 
of  an  appropriate  set  of  basis  functions: 


6(u\Cl)  —  ^  ^  tinfinju) , 
n 


(9) 


with  the  basis  set  {(/>„(•)}  selected  in  a  manner  that  allows  for  accurate  modeling  of  the 
aberration  function  [1],  When  the  basis  set  is  orthogonal  over  the  system  pupil  with 


tyn  (w)  0m  (w)  du 


L 


du 


1  n  =  m 
0  n  7^  m 


(10) 
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then 


2  _  Jy 


€  = 


[  92(u)du 

/  du 
Jy 

y'Vanom  /  (f>n(u)4>m{u)du 
„  ™  Jv 


L 


du 


=  £<4, 


(ii) 


and  the  average  square  aberration  strength  is  equal  to  the  sum  of  the  squares  of  the  coeffi¬ 
cients  {a„}. 

The  Zernike  polynomials  are  probably  the  most  commonly- used  basis  functions,  but  other 
choices  are  possible.  When  an  appropriate  statistical  model  is  available  for  the  wavefront 
errors,  for  instance,  the  Karhunen-Loeve  basis  functions  can  be  used.  Whereas  the  Karhunen- 
Loeve  expansion  has  the  advantage  of  being  optimal  in  the  sense  of  compressing  aberration 
energy  into  low-order  terms,  the  Zernike  functions  have  a  convenient  mapping  to  familiar 
aberration  modes  such  as  defocus,  tilt,  coma,  etc.  For  thin-screen  atmospheric  turbulence 
models,  however,  the  compression  properties  of  the  Zernike  polynomials  are  very  close  to 
those  of  the  optimal  Karhunen-Loeve  functions  [1].  Regardless  of  which  basis  expansion  is 
used,  the  point-spread  function  for  the  Mh  channel  can  be  modeled  as: 


hk{y,  ®)  — 


2 


exp  [j9(u-,a)]sk(y,u)du  , 


(12) 


where  the  notation  hk{y,  a)  now  shows  the  explicit  dependence  of  the  point-spread  function 
on  the  expansion  coefficients. 

In  the  remainder  of  this  section  we  consider  lower  bounds  on  the  accuracy  with  which 
the  aberration  parameters  a  can  be  estimated  from  noisy  measurements  of  the  point-spread 
functions  {hk{y\ a)}.  As  a  result  of  this  study,  insight  is  gained  into  the  following  important 
question: 


Which  modality  -  conventional  image  or  Hartmann  wavefront  sensor  image  - 
allows  for  better  accuracy  when  estimating  the  aberration  parameters? 


This  question  is  answered  through  a  detailed  analysis  of  the  Cramer-Rao  lower  bound  on 
estimation  accuracy. 
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Noise  Models  and  the  Cramer-Rao  Lower  Bound 

We  consider  a  general  scenario  in  which  light  from  a  known  point-source  is  gathered  by  an 
imaging  instrument  and  split  into  K  channels.  The  channels  might  consist  of  a  conventional 
imaging  channel,  a  Hartmann  wavefront  sensor  channel,  a  phase  diversity  channel,  or  any 
combination  of  these.  The  detection  of  optical  radiation  in  each  channel  is  fundamentally 
a  random  or  stochastic  process,  and  the  semi-classical  theory  of  photo-detection  provides 
one  means  for  modeling  this  randomness  [4].  According  to  the  theory,  the  best  possible 
detector  of  optical  radiation  will  record  photo-events  that  can  be  modeled  as  samples  from 
a  conditional  Poisson  process  with  a  rate  function  that  is  proportional  to  the  image-plane 
intensity.  Accordingly,  the  photo-event  locations  can  be  described  by  the  two-dimensional 
Poisson  counting  process  {Nk(y),  y  =  (r/1,2/2)  G  5ft2},  where  Nk(y)  denotes  the  number 
of  photo-events  that  occur  over  the  spatial  region  {y1  :  y'  <  y}  in  the  kth  frame.  The 
expectation  of  Nk(y)  conditional  on  hk(y,a )  is  then: 

E[Nk(y)\hk{y;a)}=  [  hk{y'\a)dy'.  (13) 

Jy'<y 

Detectors  of  finite  size  will  integrate  these  counts  over  spatial  regions  to  form  the  data: 

<4[p]  =  J  wkp(y)dNk(y),  (14) 

where  wkp(y)  is  a  normalized  (to  unit  integral)  weighting  function  that  models  the  spatial 
region-of- integration  for  the  pth  detector  element  in  the  kth  channel,  and  dNk(y )  can  be 
loosely  interpreted  as  the  number  of  photo-events  occurring  in  the  square  region  defined  by 
y  and  y  +  dy.  This  Riemann-Stiltjes  integral  is  interpreted  as 


/ 


wkp(y)dNk{y) 


'  0,  A4  =  0 

■A4 

Y,wkn(y>d),  A4  >  0 
,  1=1 


(15) 


where 

A4  =  /  dNk(y)  (16) 

is  the  total  number  of  photo-events  occuring  in  the  A;th  channel,  and  {yki}  are  the  spatial 
locations  for  photo-events  in  the  kth  channel.  The  conditional  expectation  of  these  data  is 

E  {dk\p]\hk(y,  a)}  =  J  wkp(y)hk(y,a)dy,  (17) 

and,  because  the  detector  regions  are  non-overlapping,  the  data  from  different  detector  el¬ 
ements  are  statistically  independent.  Because  our  goal  is  to  quantify  the  fundamental  ad- 
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vantages  of  various  sensor  modalities  for  the  estimation  of  wavefront  aberration  parameters, 
we  consider  only  the  effects  of  photon  noise  in  our  subsequent  analysis.  If  desired,  other 
noise  effects  such  as  read  noise  for  charge-coupled  device  (CCD)  cameras  can  be  included  to 
characterize  performance  limits  for  specific  cameras. 

The  collected  data  for  all  detector  elements  of  all  channels  are  denoted  by  V  =  [dk}k=l, 
and,  consistent  with  the  semiclassical  theory  for  the  detection  of  radiation,  the  conditional 
probability  mass  function  for  these  data  is  modeled  as: 


P{V;a)  =  P(dk;a) 

=  nrie*P(-Mp; «])  (Mp;  a])dfcln]  /dk[n]\, 

k= 1  P 

where 

Mp;°]  =  /  v>kv(y)hk(y\a)dy. 


(18) 


(19) 


Cramer- Rao  Lower  Bound  for  the  Estimation  of  Deterministic  Parameters 


The  Cramer-Rao  bound  is  an  information-theoretic  limit  on  the  estimation  accuracy  one  can 
expect  when  estimating  parameters  from  random  measurements  [5].  Fundamental  to  this 
bound  is  the  Fisher  information  which  quantifies  the  parameter-specific  information  con¬ 
tained  in  the  measurements.  An  important  property  of  the  Fisher  information  is  additivity 
for  independent  measurements.  That  is,  the  Fisher  information  for  the  collection  of  data  V 
is  equal  to  the  sum  of  the  Fisher  information  for  the  data  from  each  channel: 


K 


(20) 


k—  1 


where  Jk(a)  is  an  N  x  N  matrix  defined  according  to 

[d2ln  P{dk-a) 


\J  fc(fl)]nra  —  P 


dandan 


(21) 


Consistent  with  the  Poisson  photocount  model,  the  Fisher  information  for  each  channel  is 
evaluated  as: 

dhk\p\  a]  dhk[p\  a]  1 


\J l;(®)]nm  —  ^  ) 


dan  dam  hk\p\a\' 


(22) 


where 


=  _2afeIm  |  /  wkp(v)  jv  (f>n{u)Tk(u)  exp  ^j^2an^)n{u)j  sk{y,u)g*k{y,  a)dudy^  , 


(23) 
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with 

gk(y;  a)  =  Jv  Tk{u)  exp  ^  an$i>n(u)^  sk(y ,  u)du.  (24) 

Because  both  hk[p;a]  and  its  derivative  scale  linearly  with  ak,  the  Fisher  information  also 
scales  linearly  with  this  parameter. 

The  Cramer-Rao  lower  bound  makes  use  of  the  Fisher  information  to  formulate  the 
following  bound  on  the  error  covariance  for  any  unbiased  estimator  of  the  unknown  parameter 
vector  a: 

E  ^(a  —  a)(a  —  a)rj  >  [J(a)]_1,  (25) 

where  the  inequality  implies  that  the  difference  between  the  left-hand  matrix  and  the  right- 
hand  matrix  is  non-negative  definite.  As  a  consequence  of  the  Cramer-Rao  bound, 


E 


^  1  (q.jj  Gn) 
-  n 


>  trace  {[■/(a)]-1}, 


(26) 


so  that  the  Fisher  information  can  be  used  to  place  a  lower  bound  the  expected  value  of  the 
residual  aberration  strength. 


Example 

In  this  section,  we  present  an  example  in  which  performance  bounds  are  computed  for  the 
estimation  of  wavefront-error  parameters  from  a  conventional  image  and  from  a  Hartmann 
sensor.  For  this  study,  the  Hartmann  array  contained  approximately  8  subapertures  across 
the  system  pupil  as  shown  in  Figure  1.  Lower  bounds  on  the  estimation  error  for  the 
estimation  of  the  first  30  coeficients  in  the  Zernike  expansion  of  the  wavefront  were  then 
computed  for  a  conventional  imaging  system  and  a  Hartmann  sensor  system.  Because  the 
Zernike  coefficient  vector  a  is  a  random  parameter  for  turbulence-induced  aberrations,  the 
lower  bound  on  the  residual  aberration  is  also  random.  Accordingly,  we  compute  the  bound 
for  many  realizations  of  the  turbulence  parameters  for  various  ratios  of  aperture  diameter  D 
to  seeing  parameter  r0,  where  r$  quantifies  the  turbulence  strength  [1].  The  resulting  bounds 
are  averaged  over  200  realizations  and  plotted  in  Figure  2.  The  bound  values  are  normalized 
by  the  expected  number  of  detected  photons,  so  that  the  bound  for  a  particular  situation 
can  be  obtained  by  dividing  our  bound  by  the  expected  number  of  photons.  Whereas  this 
bound  is  computed  for  the  estimation  of  only  the  first  30  coefficients,  emperical  analysis  has 
shown  the  bound  to  be  independent  of  the  number  of  coefficients  that  are  actually  present 
in  the  wavefront  aberrations.  That  is,  the  estimation  for  the  first  30  coefficients  when  40 
coefficients  are  actually  present  results  in  the  same  bound  on  performance.  In  addition, 
maximum-likelihood  estimation  of  the  parameters  produces  error  variances  very  close  to 
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Figure  1:  Lenslet  array  used  to  characterize  the  performance  of  a  Hartmann  array  for  phase- 
error  estimation.  The  telescope  aperture  is  seen  to  contain  roughly  8  lenslet  apertures  across 
its  diameter. 
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Figure  2:  Cramer-Rao  lower  bounds  on  the  estimation  of  the  first  30  Zernike  coefficients 
for  turbulence-induced  wavefront  aperrations.  The  bounds  are  normalized  by  the  expected 
number  of  detected  photons. 
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these  bounds. 


Impact  on  the  Air  Force 

The  methods  and  results  presented  in  this  section  have  a  direct  impact  on  the  design  of  future 
sensing  systems  for  Air  Force  space  surveillance.  Collaborations  on  and  dissemenation  of  this 
work  has  been  facilitated  through  Bruce  Stribling  at  the  Air  Force  Maui  Optical  Station  in 
Maui. 

Dual-Channel  Signal  Recovery  from  Auto-  and  Cross-Correlation 
Data 

In  this  section  we  consider  the  inverse  problem  of  recovering  two  non-negative  signals  from 
measurements  of  their  auto-  and  cross-correlation  functions.  Problems  that  are  closely  re¬ 
lated  to  this  include  phase  retrieval,  wherein  an  object  function  must  be  recovered  from 
a  measurement  of  only  its  autocorrelation;  and  blind  deconvolution,  wherein  two  object 
functions  must  be  recovered  from  a  measurement  of  their  convolution  (or  cross-correlation). 
Whereas  these  problems  have,  in  general,  been  cosidered  separately,  common  themes  are 
seen  in  the  many  solutions  that  have  been  proposed.  In  almost  all  cases,  for  instance,  the 
solutions  are  iterative  and  require  the  utilization  of  auxilliary  information  about  the  range 
or  domain  of  the  object  function. 

The  general  problem  we  address  can  arise  in  a  variety  of  applications  ranging  from  as¬ 
tronomical  and  space-object  imaging  to  synthetic  aperture  radar  (SAR)  surveillance.  We 

propose  and  discuss  an  iterative  solution  wherein  a  sequence  of  object  estimates  is  produced 
having  the  desirable  properties  that  each  estimate  in  the  sequence  is  non-negative,  and  that  a 
measure  of  discrepancy  (I-divergence)  between  the  measured  correlations  and  the  estimated 
correlations  decreases  monotonically  with  the  number  of  iterations. 

Specifically,  if  the  unknown  signals  are  denoted  as  Oi(-)  and  o2(-),  then  the  auto-  and 
cross-correlation  functions  for  these  signals  are  denoted  by: 

rn(v)  =  5>i(z)oi(z-y),  (27) 

X 

r22(y)  =  ^o2(x)o2 (x  -  y),  (28) 

X 

and 

ru(y)  =  5Zoi(x)o2(x  -  y),  (29) 

X 

where  ra  is  the  autocorrelation  for  Oi,  r22  is  the  autocorrelation  for  o2,  ri2  is  the  crosscorre¬ 
lation  for  Oi  and  o2,  and  we  have  restricted  our  attention  to  signals  defined  on  discrete  space 
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(or  time). 


The  Signal  Recovery  Algorithm 

We  first  state  an  iterative  algorithm  for  the  recovery  of  Ox  and  02  from  the  noisy  measurements 
of  rn,  r22,  rn: 


o?ew(x)  =  o°ld(x)-^- 


E  [°?ld(^  -v)  +  °?ld(z  +  y)}  +  E  °21(V  -  y ) 


duiy) 

rniy) 


2  0\  +  O2 


(30) 


and 


or(x)  =  o^(x)^ 


where  0\  and  02  satisfy 


and 


E  Kd(z  -  *)  +  <**(* + v)]  M + E  <d(*  +  y)^ 


O !  = 


_ r22*(y) 

2O2  +  Oi 

2  E  (?/) +  E  ^i2  (^) 

y _ y _ 

20l  +  O2 

2Erf22(2/)  +  EfMy) 


(31) 


02  = 


(32) 


(33) 


202  +  0\ 

and  the  noisy  measurments  are  of  rn,  r22,  and  ri2,  are  dn,  d22,  and  d12,  respectively.  As  we 
will  soon  show,  this  algorithm  has  the  desireable  properties  that  oJew  and  c>5ew  are  always 
nonnegative  functions  (provided  that  the  initial  estimates  are  nonnegative),  and  that 


£»new  <  Dold,  (34) 

where  Daevj  and  Dold  are  the  I-divergence  discrepancy  between  the  data  and  estimates  (new 
or  old)  defined  as: 

+  E  in  +  KT(y)  -  fe(y)]| ,  (35) 
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where 


x 

To  show  that  the  I-divergence  is  non-increasing,  we  first  note  that  the  measure  D  is 
equivalent  to  the  negative  log-likelihood  function  for  situations  in  which  the  data  {dtJ(y)} 
are  independent  Poisson  variants,  with  the  expected  values  {fy(j/)}.  Therefore,  maximizing 
the  log-likelihood  for  Poisson  data  is  equivalent  to  minimizing  D.  Next,  we  show  that  the 
iterative  algorithm  is  an  instance  of  the  expectation-maximization  (EM)  method,  and  as  such 
the  sequence  of  estimates  has  the  property  of  non-decreasing  likelihood  -  or  non-increasing 

D[  6]- 

EM  Method  for  Poisson  Data 

Suppose  the  noisy  data  are  Poisson  deviates  whose  means  are  the  actual  auto-  and  cross¬ 
correlations.  In  this  case,  the  log-likelihood  is  [7]: 

L  =  ^{dn(r/)lnrn(y)  -  ru(y)} 

V 

+  ^  {d22(y)  In r22(y)  -  r22{y)} 
y 

+  5Z{di2(y)lnri2(j/)  -ri2(y)}.  (37) 

y 

Therefore,  maximizing  L  over  0\  and  o2  is  equivalent  to  minimizing  D.  To  use  the  EM 
method  for  this  process,  we  begin  by  defining  the  complete  data  as: 

dij{x,  y)  ~  Vioi^Ojix  -  y)},  (38) 

where  the  notation  is  meant  to  imply  that  dij(x,y )  is  a  Poisson  random  variable  with  mean 
Oi(x)oj(x  —  y).  The  many-to-one  mapping  from  the  complete  data  to  the  measured  or  in¬ 
complete  data  is: 

dij(y)  =  ^dij(x,y).  (39) 

X 

The  log-likelihood  for  these  complete  data  is  then: 

Led  =  Jl{dii(a;,y)ln[oi(a:)o1(a:-t/)]  -  o^o^x  -  y)} 
x,y 

+  ^22(2;,  y)  In  [o2(a;)o2(a;  —  y)]  —  o2(x)o2(x  —  y)  j 

x,y 

+  £P  12(x,y)\n[o1(x)o2{x  -  y)]  -  Oi(z)o2(:r  -  y)}  .  (40) 

x,y 
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Following  the  EM  method,  we  form  the  conditional  Expectation  of  Lc<^\ 

Led  =  53{du(ar,y)ln[oi(ar)oi(a:-y)]-Oi(a:)oi(a:-y)} 

+  E{4(*  ,  y)  In  [< o2(x)o2(x  -  y)}  -  o2(x)o2{x  -  y)} 

x,y 

+  Y2{^{x,y)\n[o1(x)o2{x-y)]-oi(x)o2{x-y)},  (41) 

x,y 


where  dij(x,y )  is  the  expectation  of  the  complete  data  conditional  on  the  measured  (or 
incomplete)  data,  and  assuming  that  the  actual  parameter  values  are  the  current  estimates 
(o°ld  and  c>2Id).  Next,  for  the  Maximization  step  we  select  Ojew  and  o!>ew  to  maximize  Lc 
From  the  Kuhn-Tucker  conditions,  we  find  the  following  conditions  for  a  maximizer: 


'(*)  = 


E  [du(x,y)  +  dn(x  +  y,y)]  +^2dn{x,y) 

y  y 


2  Eorw+E»rw 


(42) 


and 


~new 


E  [d22{x,y)  +  d22(x  +  y,y)]  +  d12(x  +  y ,  y) 

(*)  =  -* - 


“  2Eorw+E°rw 
2  2 

These  expressions  are  simplified  by  noting  that: 

dij{x,y)  =  ofd(x)of(x  - 


so  that 


E  Kd(*  -  y)  +  °fd(x + y)]  ^573  +  E  °2ld(x  - 

0?ew(x)  =  o?ld(z)- - — - EE — y - EE 

1  v  '  iv/  o  -.new  /  j\  ,  Miew/^v 

2  2^  0 1  \z)  +  2^  °2  [z) 


(43) 


(44) 


(45) 


and 


new 
°2 


(*)  =  °2ld(z) 


E  [EE 


y)  +  ofd{x  +  y) 


\*7  / 


>^w+Vr[x+v> 


r?2d(y) 


2  eeeheee) 


(46) 
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Now.  note  that  if  we  sum  these  expressions  over  x ,  we  obtain  the  follwing  equations: 


2  51  da  (y)  +  Edi2(y) 
n  _  _ v _ 

20i  +  O2 


(47) 


and 


2  dii  (y)  +  5Z  di2  (y) 


02  = 


(48) 


2O2  T  Ox 

where  Oi  =  o™™  (x)  is  independent  of  the  iteration  number.  That  is,  the  sum  of  the 
object  estimates  does  not  change  after  the  first  iteration,  and  this  sum  is  determined  by 
these  equations.  Therefore,  these  equations  can  be  solved  and  the  sums  used  for  the  resulting 
iterative  equations: 


5]  [o°ld(x  -  y)  +  o?ld(a:  +  y) 
'(*)  =  <(*)- - 


du(]ll  +  V0' 

(y)  +y 


old 


-old 

rll 


(x-y) 


dn{y)_ 

( y ) 


.old 

12 


and 


20i  +  O2 


E  I0?"1*1  -  y)  +  °f(i  +  y)]  +  E  °T‘ix  +  v) 


o5"(i)  =  o?d(x) 


-old 
r  22 


(») 


di2(y)_ 

(») 


rold 

r12 


2O2  +  Ox 


(49) 


(50) 


Because  the  EM  method  provides  a  sequence  of  estimates  that  has  the  property  of  non¬ 
decreasing  likelihood,  it  follows  that  the  resulting  algorithm  also  has  the  property  of  non¬ 
increasing  I-divergence.  Nonnegativity  of  the  estimates  follows  naturally  provides  the  initial 
object  estimates  are  nonnegative. 


Example 

Shown  in  Figure  3  are  two  object  function  and  their  auto-  and  cross-correlation  functions. 
The  object  estimates  obtained  by  using  the  iterative  method  described  in  this  report  with 
random  starting  guesses  are  shown  in  Figure  4  for  1000  iterations. 


Impact  on  the  Air  Force 

The  methods  and  results  presented  in  this  section  have  a  direct  impact  on  the  processing 
of  data  recorded  by  sensing  systems  for  Air  Force  space  surveillance.  Collaborations  on 
and  dissemenation  of  this  work  has  been  facilitated  through  David  Voelz  at  the  Air  Force 
Research  Laboratory  in  Albuquerque,  NM. 
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•’ll  r22  ri2 


Figure  3:  Two  object  functions,  and  their  auto-  and  cross-correlation  functions. 


object  1  estimate  object  2  estimate 


Figure  4:  Object  estimates  after  1000  iterations  with  random  starting  guesses. 
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